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We present a solution to the sign problem in dynamical random matrix simulations of a two- 
matrix model at nonzero chemical potential. The sign problem, caused by the complex fermion 
determinants, is solved by gathering the matrices into subsets, whose sums of determinants are real 
and positive even though their cardinality only grows linearly with the matrix size. A detailed 
proof of this positivity theorem is given for an arbitrary number of fermion flavors. We performed 
importance sampling Monte Carlo simulations to compute the chiral condensate and the quark 
number density for varying chemical potential and volume. The statistical errors on the results 
only show a mild dependence on the matrix size and chemical potential, which confirms the absence 
of sign problem in the subset method. This strongly contrasts with the exponential growth of the 
statistical error in standard reweighting methods, which was also analyzed quantitatively using the 
subset method. Finally, we show how the method elegantly resolves the Silver Blaze puzzle in the 
microscopic limit of the matrix model, where it is equivalent to QCD. 



I. INTRODUCTION 

Our knowledge about QCD at finite baryon density is 
scarce, mainly because dynamical simulations at nonzero 
chemical potential are severely hindered by the sign prob- 
lem, see Ref. pQ for a review. The problem arises because 
the fermion determinant becomes complex and can no 
longer be included in the probabilistic weight in impor- 
tance sampling methods. At small chemical potential 
the problem is mild and can be circumvented with meth- 
ods, like reweighting, Taylor expansions and analytic con- 
tinuation from imaginary to real chemical potential [JJ. 
The cost of these methods typically grows exponentially 
with the volume and the chemical potential, such that 
they quickly become unusable. For some specific models, 
tailor-made approaches were proposed which weaken or 
sometimes even solve the sign problem [2HB]. 

Because of the equivalence between QCD in the e- 
regime and the microscopic limit of chiral random matrix 
theory (RMT) [7HFJ, useful spectral information about 
the Dirac operator in QCD can be gained from RMT, 
both at zero and nonzero chemical potential. RMT 
models for QCD at finite density where proposed by 
Stephanov [TO] and Osborn [TTJ. These models not only 
allow the computation of the densities of complex eigen- 
values of the Dirac operator [TTMl3| . which were suc- 
cessfully verified in quenched lattice simulations [T4lR6| . 
but they also enable the study of the average phase of 
the fermion determinant [TTH2"T] , which characterizes the 
severity of the sign problem in dynamical simulations. 
Moreover, the relation between the strongly oscillating 
spectrum, the sign problem and the chiral condensate 
at nonzero chemical potential was investigated, and a 
new mechanism was proposed to explain chiral symme- 
try breaking at finite chemical potential, which funda- 
mentally differs from the Banks-Casher solution at zero 
chemical potential 22 , 23J . At the same time this mecha- 
nism solves the Silver Blaze puzzle [24] , as subtle cancel- 
lations ensure that thermodynamic observables are in- 
dependent of the chemical potential below the nuclear 



matter threshold. 

The aim of our work was to develop an importance 
sampling Monte Carlo method to simulate the dynami- 
cal partition function of random matrices, without using 
additional analytical input. One of the motivations is 
to learn from numerical solutions to the sign problem in 
simpler models to make progress in the solution of the 
sign problem in QCD. 

We herein give a detailed report of a subset solution 
to the sign problem in dynamical simulations of the two- 
matrix model of Osborn at nonzero chemical potential 
[TT] , which we first proposed in an earlier letter [25] . The 
principal feature of the method is the construction of 
subsets of matrices which have real and positive weights, 
but whose cardinality only grows linearly with the vol- 
ume. This positivity was stated in the initial letter and 
its proof is given in full herein. 

The positive weights allow the use of importance 
sampling to construct a Markov chain of subsets, dis- 
tributed according to the random matrix partition func- 
tion. These subset samples were used to compute observ- 
ables, like the chiral condensate and the average quark 
number density. The errors on the measurements clearly 
show that the new method is free of sign problem. This 
strongly contrasts with the results of standard reweight- 
ing methods, where the exponential blowup of the errors 
clearly signals the existence of the sign problem and leads 
to the failure of these methods. 

Interestingly, the statistical error on the average 
reweighting factors, which drives the sign problem in the 
reweighting methods, is a quantity that can also be com- 
puted using the subset method. This allowed for a quan- 
titative determination of their exponential increase and 
for a comparison of various reweighting schemes. 

We further illustrate that the subset method remains 
efficient in the parameter region where large oscillations 
in the Dirac spectrum are crucial to get the correct value 
for the chiral condensate [22j [23]. Moreover, we show 
how the subset method resolves the Silver Blaze puzzle 
in the microscopic limit of the RMT model, where it is 
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equivalent to QCD in the e-regime. 

Other approaches based on partial integrations or sum- 
mations were considered in Refs. [2"5H50] . 

This paper is composed as follows. In Sec In] we intro- 



duce the random matrix model. In Sec. Ill we present the 
subset method and give the positivity theorem, which al- 
lows for its use in Monte Carlo simulations. In Sec. lIVI we 
compare the numerical results for the chiral condensate 
and the quark number density obtained with the subset 
method and with standard reweighting methods, and ex- 
plicitly compute the reweighting factors occurring in the 
latter. In Sec. [V] we discuss some implications of the pos- 
itivity relation and show how the subset method resolves 
the Silver Blaze puzzle in the microscopic limit of RMT. 
Finally, we conclude in Sec. |VI| In the appendix we give 
a proof of the positivity theorem, details on the numeri- 
cal implementation, an overview of reweighting and some 
analytical RMT formulae. We also show how the Silver 
Blaze can be realized away from the microscopic limit. 



II. RANDOM MATRIX MODEL 

In the two-matrix model of Osborn the complex 
matrices 0i and 02 of dimension (N + v) x N are dis- 
tributed according to the unquenched partition function 

, N f 
Z= / d0id0 2 w(4> 1 )w((f>2) Y[ det £>(</>!, 2 ;^, to/), 

/=i 

(1) 

where the integration is performed over the real and 
imaginary parts of all matrix entries. The weights consist 
of a Gaussian part and a fermionic part originating from 
the Nf dynamical quarks with masses to/. Each random 
matrix has a Gaussian weight 

w^i) = ( 7 iV/^) JV ( Ar+l ')exp(-7iVtr^^), (2) 

where 7 G R + and we will adopt the convention 7=1 
to conform with the standard discussions of this random 
matrix model. An alternative choice for 7 will be briefly 
discussed in Sec. [V] Each fermion flavor contributes to 
the partition function with the determinant of its Dirac 
operator, which in this two-matrix model is given by 



D(<j>i, 2 ;^, to) = 



TO i<j) X + ^0 2 . 



for a fermion of mass to at chemical potential fi (D has 
v zero modes for m = 0, where v > without loss 
of generality). In the presence of a chemical potential 
the fermion determinant becomes complex, such that the 
fermion weights can no longer be included in the prob- 
ability distribution used in importance sampling Monte 
Carlo simulations. Standard solution methods used to 
circumvent this issue typically suffer from the sign prob- 
lem, as the work needed to make reliable measurements 
on the statistical ensemble grows exponentially with the 
volume pp. 



III. SUBSET METHOD 

Below we present a method, first introduced in 
Ref. [25] , which avoids the sign problem in dynamical 
simulations of random matrices. The main idea is to 
gather matrices into subsets which have net real and pos- 
itive weights in the partition function (JlJ. 

The subsets are constructed as follows: For any config- 
uration $ = (0i, $2) we consider a set of configurations 



fi($) = *($;#„) 



An = 0,...,N s -l\, (4) 



N 



where 9 n ) = (ipi,ip2) are orthogonal rotations of the 
seed configuration $ defined as 



ipi\ _ ( cos 9 sin6>\ f(j>i 
i>2) ~ I— sin.0 cosfly I 02 



(5) 



The initial motivation leading to the construction of 
the subsets using orthogonal rotations can be found in 
Ref. [35] ■ The product w(ipi)w(ip2) is independent of 9 
under the orthogonal rotation such that all the con- 
figurations in a subset have the same Gaussian weight, 
which we denote as W(fl). 

The random matrix partition function ([7J) can then be 
rewritten as an equivalent integral over subsets Q defined 
in Q, such that 



(6) 



Z= / dnW(n)an(n,m), 



where the fermionic subset weight, for Nf degenerate 
quarks of mass to, is given by the sum of complex de- 
terminants, 



cto(/z,to) 



n s -i 



det N fD(V n ;n,m), 



(7) 



with = 9 n ) G £1. The subset partition function 
6 ) is equivalent to the random matrix partition function 
1) as each configuration $ = (0i,0 2 ) of the random 
matrix ensemble can be used as seed of the corresponding 
subset fi( < I > ), defined in Eq. Q, and the set of all subsets 
forms an iV s -fold covering of the original random matrix 
ensemble. 

The success of the method is based on the follow- 
ing positivity theorem: For any Q constructed accord- 
ing to Q, and for arbitrary /1 < 1 and mass to, the 
fermionic subset weight an(fj,,m) defined in ([7]), i.e., the 
sum of fermionic determinants of the N s configurations 
^($; 9 n ), is real and positive if N s > NfN. 

This theorem immediately follows from the identity 



MM, to) = (1 - fl Y f(N+v/2) °n ( 0, y=f J > (8) 

which holds for any fj, and to if N s > NfN and relates 
the fermionic subset weight at chemical potential /x and 
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and effective mass = n 
the right-hand side of Eq. 



mass m to the subset weight at zero chemical potential 

y/l — [L 2 . When looking at 
_J, we observe that for zero 
chemical potential and mass £ K the Dirac matrix 
gcncrically has v real eigenvalues and N complex 
conjugate pairs to m ± zA with A € M + , such that its de- 
terminant is positive. Hence, this right-hand side is real 
and positive for fj, < 1, and so too will be the subset 
weight crfi(jj,,m) on the left-hand side. For fi > 1 the 
effective mass becomes imaginary and the subset weights 
are no longer guaranteed to be of definite sign altogether. 
In the limit m = (and v = as the determinant is zero 
elsewise) the relation becomes 



(9) 



for any [i, which is the equation given in Ref. |25j . In 
that first publication an inequality was given for the case 
m =/= 0, which is now replaced by the general identity (J8|. 

Note that for fi = 1 Eq. M yields <rn(l, 0) = 0, i.e., the 
sum of determinants exactly vanishes for /i = 1 and m — 
0. This corresponds to maximal non-hermiticity, where 
the average phase factor and the partition function are 
exactly zero, and the sign problem is maximal when using 
traditional solution methods. For nonzero mass, the limit 
of Eq. ([§ for fx -> 1 is lim^j a Q (fi, m) = N s m N f( 2N+l, '> 
and all the subsets in the partition function have identical 
fermionic weights, even though their Gaussian weights 
W{£1) will vary. 

The detailed proofs of the identities Q and ^ are 
given in Appendix [Aj 

The positivity implies that the subset weights 
W(fl)<r n (fi, m) can be used to generate subsets of random 
matrices using importance sampling methods like the 
Metropolis algorithm. As N s has to be larger than NfN 
to ensure positivity, it will be set to its optimal, i.e. small- 
est possible, value in the simulations, N s = NfN + 1. 

Using a sample of Amc subsets k = 1...Nmg, 
the expectation value of an observable O in the RMT 
ensemble is approximated by the sample average 



-, Nmc 



<>/, 



where the subset measurement is defined as 

1 N.-l 

(0) a = — V det jv ^(*„)0(* n ), 



(10) 



(11) 



n=0 



with \& n G f2 and (1)q = 1, and we omit the arguments /i 
and m from now on when no confusion is possible. The 
subset measurement is a modified subset average, which 
takes into account that the subsets are generated with a 
fermionic weight ctq, while the individual matrices ought 
to be weighed by their respective Dirac determinants. In 
the simulations reported below, the subset weights and 
subset measurements were computed exhaustively, i.e., 
the determinants of all N s configurations in the subset 
were evaluated numerically and used in Eqs. |7| and (111. 



IV. NUMERICAL RESULTS 

A. Chiral condensate and quark number density 

We applied the subset method to compute the chiral 
condensate and the quark number density in the random 
matrix ensemble, which are given by 



1 d\ogZ 



2N dm 



1 

2iV 



trZT 



and 



1 d\ogZ 
2N dfi 



\2N 



02 

A 







D~ 



(12) 



(13) 



respectively, following the derivation in Append ix |B) 

Using the block structure of the Dirac matrix (I3j) more 
efficient formulae can be derived for the numerical eval- 
uation of the fermion determinant, the chiral condensate 
and the quark number density, which are are given in 
Eqs. |B3l), (|B8| and (|B10|), respectively. 



The simulations use the Metropolis algorithm to gen- 
erate subsets fi according to their statistical weights 
PU(f2)crf2. To compute the fermionic subset weights an 
the determinants of the N s matrices in the subset are 
evaluated numerically and accumulated. It is important 
to emphasize that, even though the complex determi- 
nants fluctuate strongly when the chemical potential and 
volume increase, there is no sign problem in the com- 
putation of the individual subset weights an, as all N s 
contributing determinants are added and no statistical 
sampling is used. The subset weights are computed in a 
numericaQ but deterministic way. 

The statistical sampling comes in when the succes- 
sive subsets are generated in the Markov chain. This 
happens as follows: Start with a randomly chosen seed 
configuration and construct the corresponding set fl us- 
ing Eq. Q. Assume now that the Markov chain has 
reached a subset fit at Monte Carlo time t, then ran- 
domly choose a configuration in the current subset, gen- 
erate a new configuration by making a random step on 
each matrix entry, construct the subset corresponding to 
this new seed configuration, and apply an accept-reject 
step to the newly proposed subset to generate the subset 
fit+i. This stepping procedure is repeated until we have 
generated a large enough sample to perform the desired 
measurement. 

In our simulations, we generated iV M c = 100,000 
subsets in each Markov chain, after equilibration was 
reached. The expectation value of an observable is eval- 
uated by making a sample average of subset measure- 
ments, as prescribed in Eqs. ( 10 ) and ( 11 ). Each individ- 



ual subset measurement is computed as a deterministic 



1 This contrasts with some other partial resummation/intcgration 
methods, where partial sums/integrals in the partition function 
are evaluated analytically, and only the integral over the remain- 
ing degrees of freedom is sampled numerically. 
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FIG. 1. Chiral condensate E (top row) as a function of the chemical potential /i 2 for the subset method and the phase-quenched 
reweighting method for N — 2,4, 8. The solid line shows the exact analytical result of Eq. (D3 1. The reweighting method fails 
for ever smaller y} when N grows (some data points are negative and are left out of the semi-log plots). The corresponding 
relative statistical error e is shown in the bottom row. The error for the reweighting method grows very rapidly and should 
only be trusted as long as the method works (see top row). 



sum over N s contributions. The quantities that fluctu- 
ate statistically during the Monte Carlo sampling are the 
subset measurements (111. Successive measurements in 



the Markov chain are correlated and the number of inde- 
pendent measurements is smaller by a factor 2t, where 
r is the integrated autocorrelation time. The statisti- 
cal errors on the measurements are determined using the 
standard error formula corrected for these autocorrela- 
tions. 

To compare the subset method with standard reweight- 
ing methods, described in Appendix [Cj the simula- 
tions were repeated using quenched, phase-quenched, 
/i-quenched and sign-quenched reweighting, which are 
all expected to suffer from the sign problem ^T\. In 
rewei ghtin g methods observables are computed using the 
ratio ( C2 ) , where both numerator and denominator de- 
crease exponentially with increasing volume. The ex- 
ponential increase of the work comes from the need to 
compute these exponentially small numbers from a sta- 
tistical sampling of largely canceling contributions. The 
reweighting factors, i.e., the denominators in Eq. (C2|, 
will be discussed further in Sec. lIV Cl For the simulations 
with reweighting methods we used Nmc x N s random ma- 
trices in the Markov chains, such that the total number of 
generated matrices is the same as in the subset method. 
For the sake of clarity, we only show the results of phase- 
quenched reweighting in the figures below, as its results 
are representative for the various reweighting schemes. 

We performed simulations with one dynamical 
fermion, i.e., Nt = 1, of mass m = Q.1/2N and matrix 
sizes N = 2, . . . , 34. The mass was chosen to be small 



with respect to the magnitude of the smallest eigenvalue, 
to ensure that dynamical effects are important. 

We measured the chiral condensate given by Eq. (112 
These results were first presented in Ref. [35] • In Fig. 
the condensate is shown as a function of the chemical 
potential for matrices with sizes N = 2, 4, 8. We compare 
the results obtained using the subset method with those 
from phase-quenched reweighting and with the analytical 
results of Eq. ( D3 ) . The data, displayed in the top row, 



show that the reweighting method fails for smaller and 
smaller /i 2 as the matrix size increases, due to the sign 
problem. This strongly contrasts with the results of the 
subset method which are reliable up to much larger values 
of /i 2 and agree with the analytical predictions. 

The corresponding relative statistical errors are shown 
in the bottom row. At very small [i, when the sign prob- 
lem is not yet tangible, the error on the condensate is 
somewhat smaller for the reweighting method than for 
the subset method. This is easily explained by noting 
that at ji = the determinants are all real and posi- 
tive, such that importance sampling can be performed on 
the random matrices themselves. Sampling the partition 
function using subsets is then evidently somewhat less 
efficient. This feature persists for small, nonzero /z, but 
very quickly the exponential growth of the error in the 
reweighting method, caused by the sign problem, makes 
the method unusable. At some value of \i the error es- 
timate becomes meaningless, as the reweighting method 
completely fails. However, for the subset method, the 
relative accuracy of the measurements is nearly indepen- 
dent of the chemical potential, which confirms the ab- 
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FIG. 2. Relative error £ on the chiral condensate ver- 
sus matrix size N for various values of chemical potential 
ji 2 = 0.1,0.2,0.3,0.4,0.5. The top plot shows the results of 
the subset method, for a fixed number of subsets. The full 
curve e(N) oc y/N serves to guide the eye. As a comparison, 
the bottom plot shows the relative error for phase-quenched 
reweighting on a semi-log plot (the color coding for ii 2 is the 
same as in the top plot). 



sence of a sign problem and underscores the usefulness of 
the method. 

We also studied how the relative statistical error on the 
chiral condensate varies as a function of the matrix size 
TV for fixed values of /z 2 , and show this iV-dependence in 
Fig. |] for various values of /i 2 . For a fixed number of sub- 
sets the error in the subset method (top panel) increases 
approximately as y/~N and is independent of [i (the latter 
was already observed in Fig. [IJ. If we fix the number 
of matrices, rather than the number of subsets, the er- 
ror will increase with an additional factor y/N (as the 
subset size itself grows with N + 1), such that the over- 
all relative error will grow approximately linearly with 
N. Conversely, to achieve a constant error the number 
of subsets would have to grow proportionally to N, i.e., 
the total number of matrices should approximately grow 
as TV 2 . The bottom plot shows the same quantity for 
phase-quenched reweighting (on semi-log scale). We ob- 
serve that the error grows exponentially with N until the 
reweighting method fails and the error is no longer reli- 
able. Note that for both methods the additional cost for 
the numerical computation of the determinants is pro- 
portional to iV 3 . 

Using the same Monte Carlo algorithm we also com- 
puted the quark number density given by Eq. (13 1. The 



in Fig. [3| The data in the top row confirm that the 
reweighting method badly suffers from the sign problem 
as fi and N increase, while the subset metho d ver y nicely 
reproduces the analytical predictions of Eq. ( D4 ) . In the 



bottom row we observe that, for the reweighting method, 
the magnitude and variation of the error are very similar 
to that on the chiral condensate and clearly signify a sign 
problem. However, for the subset method, the behavior 
of the statistical error not only confirms the absence of 
a sign problem, but we observe that the relative error 
is much smaller than in the case of the chiral conden- 
sate. Moreover, even at small /x the error of the subset 
method is three orders of magnitude smaller than for the 
reweighting method, which is a very different behavior 
than for the chiral condensate. This surprising feature is 
due to the small variance of the quark number density 
over the Markov chain. This small variance can be un- 
derstood from the results derived later in Sec. |VB[ where 
we show that the subset measurement (111 of the number 



variation of the number density and its relative statistical 
error as a function of the chemical potential are shown 



density in this matrix model consists of a constant term, 
which is identical for all subsets, and a smaller term pro- 
portional to the chiral condensate; see Eq. (27 1. The first 
term does not contribute to the variance, such that the 
error on the number density is solely driven by the error 
on the latter, which explains why it is so small. From this 
we conclude that the quark number density in this model 
is especially well sampled by the importance sampling of 
subsets. The existence of subsets with large constant con- 
tribution to the quark number density is an interesting 
feature that could point to a possible search direction in 
other theories suffering from the sign problem. 

We also measured the quark number density as a func- 
tion of the matrix size N using the subset method and 
show the relative error on these measurements in Fig. [4j 
for fixed number of sampled subsets. Here again, we only 
observe a mild dependence of the error on the matrix size, 
which confirms that there is no sign problem in the subset 
method. 



B. Spectral decomposition of the chiral condensate 

The spectral decomposition of the chiral condensate at 
nonzero chemical potential is quite remarkable, and the 
ability of the subset method to cope with this intricacy 
and reproduce the correct results is yet another test for 
its viability. 

For this investigation we computed the variation of the 
chiral condensate as a function of the quark mass for 
a fixed chemical potential, both in the dynamical case 
(Nf — 1), using the subset method, and in the quenched 
case. The results are shown in Fig. [5j As was shown 
in Refs. [22, 25] the unquenched chiral condensate has a 
very different spectral decomposition depending whether 
the quark mass is outside or inside the cloud of complex 
eigenvalues of the Dirac operator, which has a width of 
about 2/j, 2 [T3]. In Fig. [6] we show typical spectral den- 
sities of the Dirac operator as the mass moves from out- 
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FIG. 3. Quark number density n (top row) as a function of the chemical potential /i 2 for the subset method and the phase- 
quenched reweighting method for TV = 2,4,8. The solid line shows the analytical result of Eq. (D4l. Again, the reweighting 
method fails for ever smaller fj? when TV grows. The corresponding relative statistical error e is shown in the bottom row. The 
error for the reweighting method grows very rapidly and should only be trusted as long as the method works (see top row). 



side to inside the eigenvalue spectrum. When the mass is 
outside this cloud, the dynamical spectrum of the Dirac 
operator is almost identical to that of the quenched case 
(see leftmost panel of Fig. [6| and the dynamical chiral 
condensate should be very close to its quenched value. 
This is confirmed in Fig. [5j where we observe that both 
curves fall together when m > 2fj?. When the quark 
mass enters the cloud of eigenvalues the quenched chi- 
ral condensate steadily drops to zero. However, as can 
be seen in Fig. [5j the dynamical chiral condensate does 
not follow this trend and remains large all the way down 
to microscopically small masses. This peculiar behavior 
was explained in Ref. [33J by the very large oscillations 
which emerge in the unquenched eigenvalue spectrum at 
z = ±m when the mass enters the cloud of eigenvalues; 
see the spectral densities plotted in Fig. [6] Subtle can- 



cellations in these large spectral fluctuations compensate 
for the decline in the quenched contribution such that 
the dynamical chiral condensate remains large when the 
quark mass is inside the cloud of eigenvalues. 

Although these cancellations were computed analyti- 
cally in Ref. [23], it seemed unlikely that the dynam- 
ical chiral condensate could be determined to a good 
accuracy through numerical simulations, at least not in 
the region where the cancellations are important, as this 
is exactly where simulations are hampered by the sign 
problem [35]. Nevertheless, Fig. [5] clearly shows that the 
subset method is able to compute the chiral condensate 
accurately, even in this critical parameter region. The 
figure illustrates how the unquenched and quenched chi- 
ral condensates move apart when the mass enters the 
cloud of eigenvalues: The dynamical condensate remains 
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FIG. 4. Relative error e on the quark number density versus 
matrix size TV for various values of chemical potential /i 2 = 
0.1, 0.2, 0.3, 0.4, 0.5 using the subset method. 



FIG. 5. Dynamical and quenched chiral condensate as a func- 
tion of the quark mass for fi 2 = 0.6 and TV = 16. The dynam- 
ical condensate was measured using the subset method. 



FIG. 6. Spectral density p(x + ly) of the Dirac operator for m = 1.2,0.4,0.1 (from left to right) with /x 2 = 0.6 and N = 16, 
computed using the formulae from Ref. [23]. For clarity, the height of the large oscillations was truncated in the figure. 



large even though the quenched value steadily decreases. 
The efficiency of the subset method remains unaffected 
when the mass enters the cloud of eigenvalues and large 
spectral fluctuations are crucial to the determination of 
the chiral condensate. A similar mechanism relating the 
spectral oscillations of the baryon number Dirac operator 
to the quark number density was recently uncovered [31] . 



C. Reweighting factors 

Another quantity that is accessible to the subset 
method is the average reweighting factor occurring in 
the denominator of the reweighting formula ( C2 ) . These 



reweighting factors are important quantities in the study 
of the sign problem in reweighting methods as they de- 
crease exponentially with increasing volume and chemical 
potential, and give rise to the exponentially growing error 
on the measurements. 

The average reweighting factor for a target ensemble 
with complex weight w simulated in an auxiliary ensem- 
ble Wq is given by the expectation value (w/wq) w , see 
Appendix [C] The direct computation of this expectation 
value in the auxiliary ensemble is obviously plagued by 
the sign problem, as it is precisely at the origin of the 
problem. However, this expectation value can be rewrit- 
ten as 



\w /, 



it'O 



J dxw(x) 
J dx wq(x) 



/wo 
\ w 



(14) 



which means that the reweighting factor in the auxil- 
iary ensemble can be computed as an inverse expectation 
value in the unquenched ensemble. Although this cannot 
be evaluated with standard methods, as the unquenched 
ensemble has complex weights, it can easily be done using 
the subset method. 

The unquenched expectation value r" 1 is then com- 
puted as the sample average ( 10 ) of subset measurements 



(111 of the inverse reweighting factor given by (assuming 
that the auxiliary weight only modifies the fermionic part 
of the action) 



\ w i n 



1 



W (*n)- 



(15) 



As the auxiliary weights Wq are real and positive, so are 
the subset measurements (15 1. Hence, the sample aver- 



ages ( 10 1 do not involve any cancellations and the average 



reweighting factor ( 14 1 can be efficiently determined by 
the subset method without encountering a sign problem, 
even though its value becomes exponentially small as the 
volume increases. 

This enables us to compute and compare the reweight- 
ing factors for various reweighting schemes, which are 
summarized in Table [I] and briefly described below: 

• quenched: the configurations are generated by di- 
rect sampling of the Gaussian weights; the average 
reweighting factor is the average fermion determi- 
nant in the quenched ensemble, 

• phase-quenched: the auxiliary ensemble is gener- 
ated using the magnitude of the complex determi- 
nants; the average reweighting factor is the average 
phase in that ensemble, 

• /i-quenched or Glasgow scheme: the auxiliary en- 
semble is generated at zero chemical potential; the 
average reweighting factor is the average ratio of 
the determinants at /i and fi = in that ensemble, 

• sign-quenched: the configurations are generated ac- 
cording to the absolute value of the real part of the 
determinant; the average reweighting factor is the 
average sign of this real part in that ensemble. This 
reweighting scheme minimizes the relative variance 
of the reweighting factors [32] ; see also [2TJ [33] . 



rew. scheme 


ferm. part of wo 


reweighting factor 


quenched 
phase-quenched 

/x-quenched 
sign-quenched 


1 

R N f 
det JV /D(0) 
\Re det N fD(fx)\ 


det N fD(n) 
exp(iNfip) 
det N fD(fj,)/det N fD(0) 
sgnRedet^D^) 



TABLE I. We list the four reweighting schemes studied in this 
paper. The auxiliary weights Wq are products of the Gaussian 
weights and a fermionic part, given in the second column. 
The corresponding reweighting factor for each configuration 
is given in the third column. The fermion determinant is 
written as det D — Re %{p . 




FIG. 7. Average two-fermion phase (e 2ie ) pq versus chemical potential in the Nf = 2 phase-quenched ens emble for m = 0.1/2JV 
and N = 2, 4, 8. The results of the subset method (blue bullets) agree with the exact result of Eq. ( |D6[ ) (solid line), while the 
direct measurement in the phase-quenched ensemble (red squares) clearly suffers from the sign problem. 
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FIG. 8. Average reweighting factor r versus chemical potential for the four reweighting schemes of Table [I] for Nf — 1, 
m = 0.1/2N and N = 2,4,8. The abbreviations in the legend are: (q) quenched, (pq) phase-quenched, (fiq) /i-quenched, (sq) 
sign-quenched. 



The aim of this measurement was to compute the 
reweighting factors for these four different schemes using 
the subset method, as described in Eqs. (10) and (151, in 



order to verify and compare their exponential decrease 
when the chemical potential and volume are increased. 

First, we verified the accuracy of the method for two- 
flavor (Nf — 2) phase-quenched reweighting, by compar- 
ing the numerical data for the average reweighting factor 
with the analytical predictions given in Appendix |D 2| 
The reweighting factor was computed using the subset 
method, as explained above, and also directly in the aux- 
iliary ensemble. The agreement between the numerical 
and analytical data is illustrated in Fig. [7] The very 
rapid decrease of the average phase factor is perfectly re- 
produced by the simulation data of the subset method, 
while the direct measurement in the phase-quenched en- 
semble fails because of the sign problem. 

In general, the average reweighting factors in the var- 
ious reweighting schemes cannot be computed analyti- 
cally. However, the subset method allows us to access 
these quantities numerically to good accuracy. We com- 
puted the average reweighting factors ( 14 ) for the four 
schemes listed in Table [I] for one flavor (Nf = 1) and 
compare the results in Fig. [8] 

These data allow us to pinpoint the onset of the sign 
problem in the phase-quenched reweighting scheme by 
first locating the /i-values where the reweighting method 
breaks down for the chiral condensate and quark num- 
ber density in Figs. [T] and [3j and then reading off the 



corresponding reweighting factors for these /i-values in 
Fig. [8] We find that the breakdown of the reweighting 
method, due to the sign problem, occurs when the aver- 
age reweighting factor drops below sa 0.01. 

When comparing the four schemes in Fig. [8j we observe 
that the Glasgow scheme has a somewhat larger reweight- 
ing factor than the other schemes. Although this could 
naively be interpreted as hinting at a weaker sign prob- 
lem, it is in reality only due to the fact that the magni- 
tude of the determinants increases with increasing /i, such 
that the values w/wq sampled in the Glasgow scheme 
are larger in magnitude than those in the other schemes. 
However, the sign problem is not actually caused by the 
size of the reweighting factor, but by its relative error, 
as it is the latter which propagates to every observable 
through Eq. (C2|. The reweighting factor itself is only 



an indicator for the sign problem, because an exponen- 
tially small value tells us that huge cancellations must 
take place. To describe the exponential problem quanti- 
tatively one has to compute the relative error e r on the 
average reweighting factor, 



2r g t 

Nmc r ' 



where 



(Rew) 



(16) 



(17) 
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FIG. 9. Standard deviation ay on the reweighting factor in various reweighting schemes as a function of the chemical potential 
for Nf = 1, m = 0.1/2N and N = 2,4, 8, computed with the subset method (blue bullets) and directly in the auxiliary ensemble 
(red squares). From top to bottom we consider the quenched, phase-quenched, /i-quenched and sign-quenched ensembles. 



is the variance of the reweighting factor in the auxiliary 
ensemble. The variance involves the second moment of 
(the real part of) the reweighting factor, which can either 
be computed directly in the auxiliary ensemble, without 
encountering a sign problem, or can be computed using 
the subset method after rewriting it as: 



Mi 



(Rew) 2 



Re w 
w 



(18) 



Rc w 



In Fig. [9] we show the standard deviation a r on the 
reweighting factor for the four reweighting schemes, 
where the second moment M2 is either computed using 
the subset method or directly in the auxiliary ensemble 
(the reweighting factor r is always computed with the 
subset method). We observe that the standard devia- 
tion is computed more accurately when Mi is directly 
calculated in the auxiliary ensemble, except for the /x- 
quenched ensemble where the subset method is much 
more accurate. Note that the latter ensemble is the only 
one where the standard deviation on the reweighting fac- 
tor substantially grows with increasing chemical potential 
and volume. 



We now merge the results of Figs. [8] and [9] to compute 
the relative standard deviation oy/r, which is propor- 
tional to the standard error e r of Eq. (fl6]). In Fig. 10 



we 

compare a r /r for the four different reweighting schemes. 
The exponential problem becomes immediately clear, as 
the explosive growth of cr r /r with increasing /1 and N 
has to be compensated in simulations by sampling an ex- 
ponentially large number of configurations to reach an 
acceptable accuracy for these reweighting schemes. The 
figure shows that, as expected, sign-quenched reweight- 
ing is doing slightly better than the other schemes [32] . 
even though it is closely followed by the phase-quenched 
scheme. For practical purposes the differences are how- 
ever not really relevant as all reweighting schemes en- 
counter the sign problem very early on. It is also in- 
teresting to note that the quenched reweighting scheme 
is not really outperformed, as we can generate uncorre- 
cted random matrices by direct sampling of the Gaus- 
sian weights (this is specific to random matrices, as the 
matrix entries are distributed independently). For the 
other schemes we use the Metropolis algorithm to gen- 
erate the configurations in the auxiliary ensemble such 
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FIG. 10. Relative standard deviation oy/r on the average reweighting factor for various reweighting schemes for Nf = 1 and 
m = 0.1/2^ versus chemical potential for N = 2, 4, 8. 



that autocorrelations have to be taken into account. The 
somewhat larger relative variance in the quenched scheme 
can then be compensated by the larger number of inde- 
pendent configurations generated for the same amount 
of work. On the other hand, the study shows that the 
Glasgow scheme performs worse than the other schemes. 

The quantitative analysis in this section confirmed that 
the sign problem is strong in all four reweighting schemes, 
and reliable numerical results can much better be ob- 
tained with the subset method, as seen in Sec. IIV A 



V. DISCUSSION 

A. Thermodynamic observables 

In this section we will see that Eq. ^ allows us to de- 
rive some interesting relations for the thermodynamical 
observables. 

For this, we first observe that the relation ^ for the 
subset weights percolates straightforwardly to the parti- 
tion function, such that 



D. Remark on the Markov chain construction 

In our simulations we computed the fermionic subset 
weights an{^,m) by evaluating and adding up the N s 
complex determinants at nonzero /i. As an alternative, 
the Markov chain of subsets can be constructed by using 
the analytic formula (|8| to compute the sampling weights 
<7n(/x,m) from crrj(0, m^). In this case one evaluates and 
sums up the fermionic determinants at [i — and effec- 
tive mass to m = to/ yl — /x 2 , and multiplies the sum with 
the suppression factor (1 — fj, 2 ^ N f( N+v ' 2 > . This sum only 
involves positive real numbers such that no cancellations 
take place. The exponential smallness of the weights at 
large fi and iV comes entirely from the suppression fac- 
tor, which is computed explicitly. As this factor is com- 
mon to all subsets, it will drop out when taking ratios 
of probabilities in the Metropolis algorithm, such that 
it plays no role when generating the relevant subsets of 
the ensemble. This implies that the relevant subsets at 
chemical potential fi and mass to are the same as those at 
fi = 0, albeit at a different, effective mass to m . Moreover, 
the relevant subsets are independent of n in the massless 
case. 

Note, that this only concerns the Markov chain con- 
struction. The subset measurements (11) still require 



the determinants at the simulated /i and to, such that the 
simulation time will increase when using this alternative 
way to construct the Markov chain, as the N s determi- 
nants have to be computed both at zero and nonzero 
chemical potential (the latter can be restricted to the 
independent subsets in the Markov chain if the autocor- 
relation times are known). 



Z Nf (fx;m) = (1 - fi 2 ) Nf(N+v/2) ZN, (0;m„) . (19) 

This relation agrees with the analytical expression for the 
partition functions derived using the method of orthogo- 
nal polynomials 23, 34 . 



Using Eq. (191 one can relate the chiral condensate at 
/i / and /i = 0. The chiral condensate for Nf = 1 is 
defined as 



1 d 

s (m;w) = — j— logZi(/x;ro), 
21\ am 



(20) 



and using Eq. ( 19 1 this can be rewritten as 
1 d 

-logZj 

i 

1 d 



2N^l - <9«V 
S(0; mjj) 



logZi (0;m^) 



(21) 



where we also used the chain rule and the definition of 
the effective mas s t% , This relation agrees with the an- 
alytical formula (D3| for the chiral condensate. 



A similar derivation can be performed for the average 
quark number density defined as (for Nf = 1) 



1 d 

n(u,:m) = —— logZi (/i;to). 



(22) 
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This can be rewritten using Eq. ( 19 1 as 

8 



n(/x; to) 



1 



1 



1 - /i 2 
_JJ_ 
1 - /i 2 
_JJ_ 
1 - /i 2 



2NJ 1 -/i 2 
v 



1 



27V 
v 

2N 
v 



log Z x (Q;m M ) 
logZi (0;m M ) 
m^E (0;m M ) 



27v 

d 

2N dm 



1+ 2iV _mS ^ ;W ) 



(23) 



where we also used the definition of to„, the chain rule, 



and Eqs. (20 1 and (21). The quark number density can 



thus be written as a sum of its massless value and a cor- 
rection term proportional to the quark condensate. This 



agrees with the analytical formula ( D4 ) for the number 
density. 

An interesting point is that the relations between the 
thermodynamic quantities at nonzero and zero chemical 
potential can also be derived at the subset level. For 
this, we note that for a thermodynamic observable (O) = 
d\og Z/dq, with q a fermionic parameter, Eq. Q yields 



(0) = 



dQW{Q) ^<r n (/i,m) 
d 

dQ.W(Q)a n (fi 7 m) — log cfq, (/j,, m) . 



(24) 



When the subsets are generated according to the sam- 
pling weights W(fl)<Jn, the subset measurements (11), 
needed to compute the sample average Eq. ( 10 ), are now 
being given by 



(0)i 



d_ 

Oq 



log <7fi(^,TO). 



(25) 



We apply this formula to compute the individual sub- 
set contributions to the chiral condensate (for Nf = 1), 
and find 



En(/x;m) 



1 d 



2N dm 
S^(0;to m ) 



log an(fi,m) 



(26) 



where the derivation is analogous to that of Eq. (21) 



Similarly we find that the contributions of the individual 
subsets to the quark number density (for Nf = 1) are 
given by, 



na(n;m) 



J^_d_ 
2Ndfi 
A* 



l-/i 



log (7n(n,m) 

1 + z^tz - toE (/x; m ) 



v 

2N 



(27) 



where we followed the same steps as in Eq. (23). The 



subset relation (27) is important in the analysis of the 
statistical error in Fig. [3] All subsets give a large com- 
mon contribution to the number density, which leads to 



a small relative variance of the quark number density. 
Moreover, for m = all subsets give the same contribu- 
tion to n, and the error on the measurement vanishes for 
the subset method, i.e., in the massless case a single sub- 
set would suffice to compute the correct number density. 
Note that the large constant contribution to the quark 
number density is a non-trivial feature of the subset con- 
struction, which cannot be identified in the contributions 
of the individual random matrices. 



B. Microscopic limit and Silver Blaze puzzle 

The equations derived above also show how the Silver 
Blaze puzzl^] is resolved in the subset method. To be 
equivalent to QCD the microscopic limit of RMT has to 
be considered, where rh = 2Nm and jl 2 — 2Nfi 2 are kept 
fixed when N oo. In this limit Eqs. (21 ) and (23) lead 
to 



m) = £ (0; rh) and n(/t;m)=0, (28) 
where we introduced the microscopic limits 



/(/t;m)= lim f(p,/V2N;rh/2N). 



N 



(29) 



The chiral condensate and the quark number density are 
thus independent of the chemical potential in the micro- 
scopic limit of RMT. Even more, following Eqs. (26) and 



( 27 1 we observe that the contribution from each individ- 



ual subset to the thermodynamic quantities is indepen- 
dent of n in the microscopic limit, such that the Silver 
Blaze puzzle is in fact already resolved at the subset level. 

Note that the prefactor in Q generates an exponen- 
tial factor in the microscopic limit of the RMT partition 
function, as 



lim (1 - ii 2 ) N ' N = lim 1 



exp 



2N 

NjfL 2 



N f N 



(30) 



Interestingly, this exponential factor is already generated 
within each subset individually as it originates from (JsJ) . 
In the RMT model the free energy density is defined as 



F(fi\ m) = - — log Z (fj,; to) 



(31) 



which also becomes independent of fj, in the microscopic 
limit because 

P(fi; to) = lim + F(0; to) = F(0; to), (32) 

N^oo 4iV 



2 The Silver Blaze puzzle in QCD refers to the fact that, at zero 
temperature and for a chemical potential less than approximately 
one third of the nucleon mass, the free energy and the thcrmo- 
dynamical observables are independent of /i 1241 . 
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FIG. 11. Chiral condensate E versus microscopic chemical po- 
tential p, 2 for rh — 0.1 and increasing values of TV. As N gets 
larger the numerical results converge towards the analytical 
microscopic value of Eq. ( D5 1 . 



where we used Eqs. (19) and (30 1. Although the factor 



( 30 1 does not occur in the partition function of chiral 



perturbation theory and is an artifact of the random ma- 
trix model, it is not relevant when discussing its universal 
properties, as it leaves the microscopic eigenvalue corre- 
lations unchanged [3TJ [35] . 

The subset method also allows a numerical investiga- 
tion of the convergence towards the microscopic limit, 
as N — > co with fixed microscopic parameters p, and rh. 
This is illustrated in Fig. [IT] where we plot the chiral 
condensate as a function of the microscopic chemical po- 
tential ft 2 for different values of N . As N increases the 
chiral condensate converges to its microscopic limit given 
by Eq. (p5|). 



It is interesting to note that choosing 7 = 1 — ri 2 (as- 
suming /1 < 1) from the onset in the Gaussian weights |2| 
not only gets rid of the spurious factor (30), but makes 



the partition function independent of the chemical po- 
tential. This was proven using the method of orthogo- 
nal polynomials in Ref. [UJ, but it also directly follows 
from the subset relation {[sj) as is shown in Appendix [E] 
This modification has the salient feature that the Sil- 
ver Blaze is then satisfied for any TV, even away from 
the microscopic limit. The /i-independence of Z does 
not, however, alleviate the sign problem of the model, as 
the fermion determinants still exhibit huge fluctuations 
when the quark mass enters the cloud of eigenvalues. In 
fact, we know that these large oscillations are essential to 
resolve the Silver Blaze puzzle at large chemical poten- 
tial. The efficacy of the subset method to solve the sign 
problem remains intact as this choice of 7 merely cancels 
the prefactor in ( 19 ) and rescales the fermion mass (see 
Eq. (E3|). The results for this 7 can easily be related 
to those previously computed with 7 = 1, and therefore 
we will not present explicit numerical results for this al- 
ternative choice. One observation is that the constant 



contribution to the quark number density (27), which 



becomes large far from the microscopic limit, will cancel 
out in this case. However, even though their results dif- 
fer at finite N, both values of 7 yield the same universal 
limit. 



VI. CONCLUSIONS 

In this paper we have presented a solution to the 
sign problem in dynamical simulations of the two-matrix 
model of Osborn at finite chemical potential. The ran- 
dom matrices are gathered into subsets, which have real 
and positive fermionic weights while their cardinality 
only grows linearly with the matrix size. A detailed proof 
of the positivity theorem for the subset weights was given. 

The positive subset weights make it possible to sam- 
ple the partition function with an importance sampling 
Monte Carlo method and generate a Markov chain of rel- 
evant subsets. As the chemical potential and the matrix 
size increase, the weights of the subsets rapidly decrease, 
but without causing a sign problem in the simulations. In 
contrast to standard reweighting methods the large can- 
cellations, inherent to simulations at real chemical po- 
tential, do not happen through statistical sampling of 
the ensemble but are confined inside the subsets, where 
the weights and measurements are computed in a deter- 
ministic way from a small number of contributions. The 
ensuing subset measurements, which are used to compute 
the sample averages, do not suffer from large statistical 
fluctuations so that the standard error on the simulation 
result is well under control. 

The method was used to compute the chiral condensate 
and quark number density accurately over a large param- 
eter range, showing that the subset method has no sign 
problem, even in regions where the reweighting meth- 
ods are unusable. The method is especially well suited 
to compute the quark number density in this model, as 
the statistical error on this quantity is extremely small, 
which was understood from analytical considerations. 

The subset method also enabled us to compute the 
reweighting factors, appearing in the standard reweight- 
ing methods, and their relative standard errors. This ex- 
plicitly revealed the exponential increase in the work re- 
quired by these reweighting methods, signaling the pres- 
ence of the sign problem. 

We also showed how the positivity relation resolves the 
Silver Blaze puzzle in the microscopic limit of the random 
matrix model, where it is equivalent with QCD, and how 
this mechanism already works at the subset level. 

The important question whether the subset method 
can be applied to physical systems and ultimately to 
QCD itself has not yet been answered and is the focus of 
current research. 
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Appendix A: Proof of the positivity theorem 

In this appendix we prove the positivity theorem for 
the subset weights, which is at the basis of the subset 
method. We first prove the massless relation Q for a 
single fermion, before extending it to an arbitrary num- 
ber of Nf massless fermions. Consequently we generalize 
the identity to the massive case and prove relation (|8|), 
first for one massive fermion and finally for Nf fermions 
with degenerate mass m. 



1. One massless fermion 

Using its block structure, the determinant of the Dirac 
matrix ([3| for a massless fermion, with v = or neglect- 
ing the zero modes, is given by (see Appendix [B]) 



det D e (fi) = det Qe(fi), 



(Al) 



where we introduce, for brevity, Dg = D 9)) for 
a rotated configuration ^(&;9) — (tpi,ip2) defined in 
Eq. |5|, and Qg = —BA is an N x N matrix with 



lip! + flips 
■ I1p\ + fllpl 



(A2) 



When expanding the product in Qg we find 

Qg(fi) = ipltpi - fi 2 ip\ip2 ~ + Wl) ' ( A3 ) 

and after substituting the definition ^ this becomes 

Qg (ft) = (cos 2 9 a + sin 2 9 b + sin 9 cos 9 c) 

- ft 2 (sin 2 9 a + cos 2 9 b - sin 9 cos 9 c) 



ifi((c 



sin 



(A4) 

2 sin # cos # (b — a)) , 



with the N x N matrices a 

AT, 



2 4>2 aud c 



0102 + 0201- After gathering the terms in a, b and c this 
can be rewritten as 



?(/•*) = f a a + fob + f c c 



with 



f a = cos 2 9 — jj 2 sin 2 9 + lifi cos 9 sin 9 
fb = sin 2 9 — fj 2 cos 2 9 — 2ifi cos 9 sin 8 
f c = (1 + /i 2 ) sin 9 cos 9 - ifi(cos 2 9 - sin 2 9) 
which, interestingly, can be further simplified to 
' fa = (cos 6 + ifi sin 0) 2 
' fb = (sin# — ifi cos 9) 2 

fc = (cos 9 + ifi sin 9) (sin 9 — ifi cos #) 



(A5) 



(A6) 



(A7) 



Using the Leibniz formula for determinants, Eq. (Al) 
can be written as 



det Dg(fi) 



N 

E 



»1)»3j— >1JV 



Qi 4l Q 



2/v 



(A8) 



where , 



is the antisymmetric Levi-Civita symbol 



and are the entries of Qg. Each term in the sum 
is a product of N matrix components, which, according 
to Eq. ( |A5[ ), can be written as Qij = f a aij + fb bjj + 
f c Cij. The coefficients f a , f b , and f c , given in §A7\, are 



functions of /i and 9, which are independent of the indices 
i and j. After expanding all the products in (A8) and 



gathering terms with equal powers of f a , fb and f c , the 
determinant can be written as 



det Dg(fi) 



N 
i,3=0 



H${a,b,c), (A9) 



where 



t»An;9) = flflf c 



i fj fN—i—j 



(A10) 



and Hfj (a, b, c) is a sum of signed products, each con- 
taining i, j and N — i — j components of the matrices a, 
b and c, respectively, which is implicitly defined by iden- 



tifying Eqs. (A8) and (A9). The /i and 9 dependence of 



the determinant ( A9 ) is completely contained in the co- 



efficients tE defined in Eq. (|A10b. Using Eq. (|A7|) these 



coefficients can be simplified to 

1%(jj.;9) = (sin0 - ifi cos 9) a (cos 9 + ifism9) p , (All) 

where a = N — i + j and j3 = N + i — j , and thus 
a + ft = 2N with < a, (3 < 2N. 

To prove the conjecture Q it is sufficient to prove that 
the identity holds term by term in the formula (A9), i.e, 



(iVfE^M, (A12) 
n 



for all < i,j < N, where we denoted the subset sum- 
mation as 



AT.-1 



n 



f(9) = fi^/Ns), 



(A13) 



with the subset defined in Eq. Q. 

To prove Eq. (A12) we first convert the trigonometric 
functions in Eq. (All I into complex exponentials using 



cos f = 



which leads to 



e tf> + e 



smf = 



[(1 



iO 



hM)e' 
+ {l-fi)e 



2i 

(l-n)e-* e ] 

910 



(A14) 



(A15) 



We apply the binomial formula to expand both powers 
and find 



%2N 



X(l 



a 13 

EE( 

fc=0 £=0 

M) 



Q 



ZN-k-t e 2i(k+l-N)d 



• fl) k+t 

(A16) 
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where we also used that a + (3 = 2N. The 9 dependence 



of Eq. (A16) is completely contained in the exponential 



function and its subset sum (A13) is 



N s -1 



(A17) 



n=0 



where q = k + £ — iV is an integer g [-N, N]. Sq is a 
sum over the q-th powers of all the A^-th roots of unity. 
It can be computed by writing it as a geometric series, 



N s -1 



(A18) 



n=0 



We distinguish two cases, depending whether q/N s is in- 
teger or not. For q/N s ^ Z we have e 2m q/N s ^ anc j 
the sum of the geometric series is 



Sn = 



1 -e 



2-Kl -j 



= 0, 



(A19) 



Nf massless fermions 



To extend the theorem to Nf > 1, the determinant 



(A9) has to be multiplied Nf times so that many more 
terms are generated. Nevertheless, after exponentiation 
the global structure of the fermionic weight remains simi- 
lar to that of Eq. ( A9 ) and can be written as the following 
linear combination 

NfN 

det N W e ((j,) = £ t?/ N (n;9)H?/ N (a,b,c), (A25) 



.NtN 



where t i / is defined in Eq. (AlOl and the implicitly 



defined function H only depends on the components of 
the matrices a, b and c, which were defined right after 
Eq. ( A4 ) . The remainder of the proof is identical to that 
for Nf — 1 with N replaced by NfN, which eventually 
leads to 

= (1 - M 2 )^ (A26) 



which vanishes as q is integer by definition. If q/N s € Z 
then e 27rt i/ N = — i and the sum (A18) can be computed 
explicitly, 



Sr. 



N. 



(A20) 



Together with Eq. (A25) this proves the conjecture ^ 
for arbitrary Nf\ 

£det^A,( M ) = (1 - f, 2 ) N f N J2^t Nf D e (0). (A27) 



If we take N s > N, the ratio q/N s will be non- integer 
and Sq zero for all q, except for q = 0, as q is an integer 
varying from — N to N. Hence, when summing Eq. ( A16 ) 
over the subset, only terms for which k + £ — N will 
survive. After considering both cases a. ^ j3, the subset 
sum can be written as 



(A21) 



where 



uj, 



2 2N 



N~ 



(zsgnA)^£(- 



fc=0 



N~ 



N+ 
N -k 



(A22) 



with N ± = N ± | A| and A = i - j. 

As uJij is independent of fj,, Eq. ( A21 1 immediately im- 
plies that 



^if(^) = (l-^)^tf(0; 



(A23) 



Because this relation holds for any i and j in Eq. (A9), 



it also holds for the sum over these indices, which proves 
the conjecture Q for Nf = 1, i.e., 

£ det De(fi) = (1 - H 2 ) N £ det D e (0). (A24) 



One massive fermion 



For the massive case we use the determinant block for- 



mula (B3) 



det Dg([i, to) ~ m v det Qe(n, m), 



(A28) 



with Q e = m 2 - BA, and A and B defined in Eq. ([A2 



To prove the conjecture (|8| for Nf — 1 we first study 
detQg- If we repeat the steps leading to Eq. (A5| we 
now find 



QeO, to) = to 2 + f a a + f b b + f c c, 



(A29) 



where / , fb and f c are defined in (|A7j). After using the 
Leibniz formula ( A8 ) and expanding all the products, as 



we did before when deriving Eq. ( A9 1 for the massless 
case, we now find 



N 



det Q (ji,m) = £ t^ h {^ 9) m 2h H^( a , b, c), 

h,i,j=Q 

(A30) 



where t-- is defined in (A10) and contains the full 9 
l j 1 



and [i dependence, while H is implicitly defined by ( A30 ) 



and depends on the components of the matrices a, b and 
c. We can repeat the whole argument of Sec. A 1 on 
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the subset sum of the ^-coefficients, after replacing N by 
N — h. In analogy to ( A23[ ) this leads to 



= (l-MT^E^M)' (A31) 
n n 

This identity still depends on the summation index h, and 



after substitution in the subset sum of ( A30 ) we find 
^det Qg(/j,,m) 



N 



Q h,i,j=0 



N 



(O;0) 



(i-mTE EC 

O h,i,j=0 



(l-MTE det ^ (°>^= 
a V v 1 



,2/j 



jV-/i 



(a,6,c) 



(A32) 



where the last equation is easily derived by inspection, 



after setting /i = and replacing to by m/yl — [j? in 
Eq. (|A30|. After multiplying Eq. (|A32| with to" and 



using Eq. ( A28 ) , we find 



det -Dg(/i, 

= (l- Ai 2 ) JV TO^detQ, fo,-^= 

= (i- M T +!//2 E det ^ (°'-7==i 



(A33) 



This proves the conjecture Q for JV/ = 1: In the massive 
case the fermionic subset weight at chemical potential fi 
and mass to is related to the weight at /i = and effective 
mass m/y/l — fi 2 . Note that for /i > 1 the effective mass 
on the right-hand side becomes imaginary, and the sub- 
set weights, although still real, can be either positive or 
negative, without having a definite sign. Therefore, im- 
portance sampling of subsets can only be used for ji < 1, 
which is the relevant region when relating the random 
matrix model to QCD. 



4. Nf massive fermions 

The proof in the previous section can be generalized to 
an arbitrary number of degenerate flavors, in exactly the 
same way as was done in Sec. |A 2| for the massless case. 
Equation ( |A32[ ) is then replaced by 



Edet^Q e ( M ,m) 



(i- Ai 2 )^E detAr/ ^(°'y ! ^)' ( A34 ) 



and after multiplying with m NfU and using Eq. (A28l, 
we finally find 



Edet^A)(/i,m) (A35) 

Q 

a \ V 1 -^ 2 / 

which proves the conjecture ^ for arbitrary Nf. 

Appendix B: Implementation 

Below we describe the numerical implementation of the 
computation of the determinant, chiral condensate and 
quark number density. We write the Dirac matrix ^ as 



D = 




(Bl) 



where A — i<j)\ + fxfa is an (N + v) x N matrix, B = 
i<fi\ + i.uf)\ is N x (N + v) and to at = diag(m . . . m) is a 
diagonal N x N matrix. 

a. Determinant. It is well-known that the determi- 
nant of such a block matrix can be computed as 

det D = det [tojv +1/ ] det [tojv — Bmjf\ V A\ (B2) 

when m 0. Note that this product of determinants 
cannot be merged as their arguments have different di- 
mensions. This can be simplified to 



det.0 = to" detQ, 



where we defined 



Q 



N 



BA. 



(B3) 



(B4) 



Note that the factor to" is reminiscent of the v exact zero 
modes of the massless Dirac matrix. In the massless case 
the determinant simplifies to 



detD = det[-BA], 



(B5) 



for v = (or for v ^ if we neglect the zero modes). Nu- 
merically, det Q was computed using an LU-factorization, 
as this is more efficient and accurate than using a full di- 
agonalization. 

b. Chiral condensate. The chiral condensate is given 
by (for Nf — 1) 



S = 



1 dlogZ 
2N dm 

1 1 
2NZ 
1 1 
2NZ 
1 



1 1 f ddetD 
2NZ dm 
3D 



det D tr 

dm 

detD trD" 1 



D~ 



2N 



trZT 



(B6) 
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where we denoted f,=J g?0ig?0 2 w(0i) w((f>2)- The in- 
verse of D can be computed efficiently using its block 
structure, which yields 



D~ 



-Q- X B 



-AQ- 



rni 



(B7) 



This formula is easily verified as DD 



1. When tak- 



ing the trace of Eq. (B7) to compute (B6) this can be 
further simplified as 



1 

2N 
1 



tr 



1 



(1 



N- 



■AQ~ l B) 



tr [mQ- 1 ] 



2Nm 
1 



2Nm 
v m 



2Nm N 



(N + v + tr [(to 2 + BA)Q- 1 }) 
(N + v + ii [(2m 2 - Q)Q- 1 ]) 
(trQ- 1 ), 



(B8) 



where we also used the definition (B4) of Q. The term 



v/2Nm originates from the zero modes of the massless 
operator and is sometimes omitted. The inverse Q^ 1 can 
be computed using the LU-factorization of Q, which is 
already available from the evaluation of det Q. 

c. Quark number density. The average quark num- 
ber density is given by (for Nf = 1) 



1 d\ogZ 



ddct.D 




\2N 



(B9) 



Using the block-inverse (B7| this simplifies to 



1 

'2N 



tr 



b\A + B(j) 2 )Q- 



(B10) 



Appendix C: Reweighting 

Reweighting methods can be used to perform Monte 
Carlo simulations when the weights are complex and the 
ensemble cannot be directly sampled with importance 
sampling. 

The ensemble average of an observable y(x) in an en- 
semble with weights w(x) is defined by 



J dx w{x)y(x) 
J dx w(x) 



(CI) 



In the reweighting method one introduces an auxiliary 
ensemble with weights wq(x) and rewrites the previous 
equation as 



(C2) 



Jdx w (x)^- )y (x) _ \w y 



If the weights wq are chosen to be real and positive, 
the auxiliary ensemble can be sampled using importance 
sampling methods and the ratio of expectation values 
in Eq. (C2) can be evaluated in a Monte Carlo simula- 
tion. Typical examples for wq are the quenched, phase- 
quenched, ^-quenched and sign-quenched ensembles. 

Reweighting methods typically suffer from an overlap 
problem, when the relevant configurations in the target 
and auxiliary ensembles do not coincide. More impor- 
tantly, when the target weight is non-positive one en- 
counters the sign problem, as the work needed to make 
reliable measurements on the statistical ensemble grows 
exponentially with volume and chemical potential be- 
cause it involves the computation of exponentially small 
reweighting factors (w/wq) Wq from a statistical sampling 
of largely canceling contributions [1] . 



Appendix D: Some analytical results 

In order to verify some of our simulation data, we quote 
a couple of known analytical results. 



1. Nf = 1 observables 

The Nf = 1 partition function Z\ can be expressed in 
terms of orthogonal polynomials as [53J 153] 



Zo 



N 



N\L V N { - 



Nm 2 



(Dl) 



where L V N are generalized Laguerre polynomials of order 
v and degree N. Their derivatives are given by 



-L v N {z) 



dz~" x ~' 

such that the chiral condensate (p0| is given by 



E(/z;m) 



-u+l ( Nm 2 \ 
J N-1 \ 



2Nm 



1 - U 2 tv { Nm 2 



(D2) 



(D3) 



The quark number density ( 22 ) can be computed analo- 
gously, yielding 



n(/z; m) 



l- f i 2 



9 J V + 1 



Ni. 



(D4) 
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In the microscopic limit, where rh = 2Nm and fi 2 — 
27V /x 2 are kept fixed while taking N —> oo the chiral 



condensate (D3) and quark number density (D4) become 



* - I' (m) 

E(/t; rh) — v - and u(/t;m) = 0, (D5) 
I v (m) 

where I v is a modified Bessel function. 

2. Phase-quenched reweighting factor for Nf = 2 

The reweighting factor, i.e. the denominator in 
the reweighting formula (C2), in the phase-quenched 



reweighting scheme is nothing but the average phase of 
the fermion determinant in the phase-quenched ensem- 
ble and can be computed analytically for Nf = 2. This 
average phase can be written as [15] 



2i9\ 



(det 2 D) Nf=0 
(\detD\ 2 ) 



N f =0 



(D6) 



Both, numerator and denominator are quenched expec- 
tation values of products of characteristic polynomials 
and their complex conjugates, which can be computed 
analytically using the method of orthogonal polynomials 
[34j . For random matrices of size N and topology v this 
yields, 



(det 2 D)jv /=0 



1 

2m 



det 



p u N (m;ii) p u N+1 (m;n) \ 

\d m p N {m;^) d m p u N+1 {m;^)J 

(D7) 



and 



(\ d etDf) N ^ = r M t mm; " )l2 



k=0 



-*(A*) 



(D8) 



with orthogonal polynomials 
,2\ k 



N 



k\Lt 



1-f 



(D9) 



and normalization factors 



Substitution of Eqs. rtD7h and (D8| in Eq. (D6) yields 



the average phase in the phase-quenched ensemble 



Appendix E: Silver Blaze at finite N 

Below we show that the partition function ([!]) can be 
made independent of /i with a judicious choice of 7 in the 
Gaussian weights pi). For arbitrary 7 e R + we rescale 
the matrices in using §\ = y/jefri, such that 



Z = C 



^exp[-JV(tr0iVi+tr^V 2 )] 



Nf 

n detD ( 
/=i 

2N(N+v) 



W7,<t>yV7^,™ f ), (El) 



and we also took into account 



where C = (N/ir) 
the Jacobian of the transformation. The 7-dependcncc 
has thus been shifted from the Gaussian weights to the 
Dirac matrix, and the partition function now looks like 
a conventional 7 = 1 partition function, albeit with a 
modified Dirac matrix. From the structure of the Dirac 
matrix (J3| we see that the scaling of the fields can be 
shifted to the mass, as 

(E2) 

Recalling that the Dirac matrix has dimension 27V + v, 
the partition function becomes, 



Z =C 



b' 2 exp[-iV(tr^Vi+tr^V2)] 



X 7 



-N f (2N+v)/2 



YldetDM.fefryfymf), (E3) 



We now look at the subset sum ^ of fermion determi- 
nants for Nf degenerate quarks in the partition function 
( E3 ) . Using the relation d8| we find that these fermionic 



subset weights satisfy 

■ _ 2X N f (N+u/2) 
1 



<Jn 



(E4) 



Clearly, if we choose 



1 — 11 2 (assuming fx < 1) 



(D10) the right-hand side of (E4) is independent of /i. From 



Eq. ( E3 1 we can then immediately conclude that Eq. ( 19 ) 



is now replaced by 

Z Nf m) = Z Nf (0; m) , 



(E5) 



such that the Silver Blaze is not only satisfied in the 
microscopic limit, but also for any finite 7V away from it. 
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